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ABSTRACT 


The  Complex  Ambiguity  Function  (CAF)  allows  simultaneous  estimates  of  the  Time 
Difference  of  Arrival  (TDOA)  and  Frequency  Difference  of  Arrival  (FDOA)  for  two 
received  signals.  The  Complex  Ambiguity  Function  Geo-Mapping  (CAFMAP)  algorithm 
then  directly  maps  the  CAF  to  geographic  coordinates  to  provide  a  direct  estimation  of  the 
emitter’s  position.  The  CAFMAP  can  only  use  short-duration  CAFs,  however,  because  the 
collector  motion  changes  the  system  geometry  over  time.  In  an  attempt  to  mitigate  this 
shortfall  and  improve  geolocation  accuracy,  the  CAFMAP  takes  multiple  CAF  snapshots 
and  sums  the  amplitudes  of  each.  Unfortunately,  this  method  does  not  provide  the  expected 
accuracy  improvement,  and  a  new  method  is  sought. 

This  thesis  reformulates  the  equations  used  in  computing  the  CAF,  in  order  to  ac¬ 
count  for  the  collector’s  motion,  and  uses  the  results  to  derive  a  new  CAFMAP  algorithm. 
This  new  algorithm  is  implemented  in  MATLAB,  and  its  results  and  characteristics 
analyzed.  The  conclusions  are  as  follows:  the  new  algorithm  functions  as  intended, 
removes  the  accuracy  limitations  of  the  original  method,  and  merits  further  investigation. 
Immediate  future  work  should  focus  on  ways  to  reduce  its  computation  time  and  modifying 
the  algorithm  to  account  for  3-Dimensional  reality,  non-linear  motion  of  collectors,  and 
motion  of  the  emitter. 
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Executive  Summary 


Knowing  the  precise  geographical  position  of  adversaries  and  their  equipment  is  critical  to 
our  ability  to  take  action  against  them.  Therefore,  the  accuracy  of  Department  of  Defense 
(DOD)  Signals  Intelligence  (SIGINT)  systems  is  of  the  utmost  importance.  These  systems 
make  use  of  a  variety  of  methods  to  perform  geolocation,  each  with  their  strengths  and 
weaknesses  [1],  [2],  [3].  This  thesis  is  an  attempt  to  improve  upon  the  accuracy  of  one  of 
these  methods,  called  the  Cross  Ambiguity  Function  Geo-Mapping  (CAFMAP)  [4]. 

The  CAFMAP  method  may  be  used  when  two  spatially  separate,  moving  collectors 
gather  independent  copies  of  the  same  emitted  Radio  Frequency  (RF)  signal.  The  Com¬ 
plex  Ambiguity  Function  (CAF)  of  these  two  signals  is  then  calculated,  which  generates 
a  surface  in  the  Time  Difference  of  Arrival  (TDOA)  and  Frequency  Difference  of  Ar¬ 
rival  (FDOA)  domains.  This  surface  peaks  at  the  actual  TDOA  and  FDOA  values  em¬ 
bedded  in  the  signals  by  the  physical  geometry  of  the  collector-emitter  system.  In  the 
CAFMAP  method,  this  surface  is  mapped  to  geographic  coordinates  and  provides  a  three- 
dimensional  (3D)  map,  peaking  at  the  estimated  location  of  the  RF  emitter.  As  shown  in 
Figure  1,  this  method  is  effective,  and  is  particularly  good  at  locating  multiple,  co-channel 
emitters  in  close  proximity  [4] . 

X  10* 
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Figure  1.  Results  of  Hartwell’s  CAFMAP  algorithm. 

The  drawback  to  the  CAFMAP  method,  however,  is  that  it  cannot  be  used  with  signals 
collected  over  long  intervals.  This  is  due  to  the  changing  collector-emitter  geometry  over 
time,  which  causes  a  smearing  effect  in  the  CAF.  When  mapped  to  geographic  coordi¬ 
nates,  this  smearing  manifests  itself  as  noise  in  the  form  of  spurious  peaks  in  the  map. 
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These  spurius  peaks  increase  the  map’s  variance  and  confuse  the  geolocation  estimate.  An 
example  of  this  is  shown  in  Figure  2.  As  the  signal  lengths  increase,  the  peaks  become 
larger  and  the  noise  floor  smaller,  but  when  the  signals  become  longer  than  ~  2^^  samples, 
the  noise  becomes  quite  noticeable. 

This  work  reformulates  the  CAF  equations  used  in  the  CAFMAP  algorithm,  such  that 
the  collector-emitter  geometry  is  updated  with  each  successive  sample  of  the  collected  sig¬ 
nals.  This  allows  for  much  longer  signal  collection  and  integration  times  without  smearing, 
which,  in  turn,  improves  the  geolocation  accuracy  of  the  CAFMAP  method. 

The  CAF  is  mathematically  defined  as 

A{tJ)=  I  +  (1) 

Jo 

where  signal  collection  starts  at  time  t  =  0,  and  ends  at  time  t  =  T.  Furthermore,  when 
applied  specifically  to  geolocating  a  received  signal: 

5i(t)  =  analytic  signal  received  at  collector  1, 

S2{t)  =  analytic  signal  received  at  collector  2, 

T  =  time  difference,  and 
/  =  frequency  offset.  [5] 

The  magnitude  of  the  CAF  is  then  maximized  when  T  and  /  are  equal  to  the  TDOA  and 
FDOA  values  inherent  in  the  two  received  signals.  These  TDOA  and  FDOA  values  are 
defined  as 


T  = 

l^2|-|0| 

and. 

(2) 

c 

k 

’vc2  ■  r2 

Vci  ■  n ' 

(3) 

c 

.  Inl 

|o| 

where  fc  is  the  transmission  carrier  frequency,  c  is  the  speed  of  light,  and  the  other  factors 
are  illustrated  in  Figure  3  [1],  [2]. 

By  allowing  the  f  and  v  factors  in  (2)  and  (3)  to  be  time  varying,  changing  by  the  known 
collector  parameters,  and  subsituting  the  result  back  into  (1),  we  arrive  at  a  new  version 
of  the  CAF,  that  now  accounts  for  the  changing  geometry  over  time.  The  result,  when 
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Figure  2.  Hartwell  CAFMAP  results  for  varying  signal  lengths. 
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Figure  3.  Example  geolocation  geometry. 
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converted  to  the  digital  domain  for  implementation,  is  shown  in  (4).  The  implementation 
of  this  modified  CAP  equation  in  the  CAFMAP  method  is 


A{x,y)  = 


no+fsT-l 

n=n() 


n-fROUND 

e 

Vc2-f2(nTs)  Vq-rpnTi) 

|r2(«7i)|  \f\{nTs)\ 

V  cTs  )_ 

(4) 


and  is  herein  referred  to  as  the  Moss  algorithm.  In  (4),  the  faetor  j  must  be 

rounded  to  the  nearest  integer  because  it  is  in  the  digital  domain,  and  corresponds  to  an 
offset  in  the  sample  index.  As  such,  fractional  offsets  are  not  allowed  as  they  would  result 
in  programming  errors. 

The  results  of  the  Moss  algorithm,  for  a  simple  ease,  are  shown  in  Figure  4.  For  eom- 
parison,  the  results  of  the  original  Hartwell  algorithm,  for  an  identical  scenario,  are  shown 
alongside.  It  is  seen  that  the  Moss  algorithm  is  still  able  to  geolocate  the  emitter.  As 
an  additional  benefit,  it  demonstrates  a  smoother,  higher  resolution  map  than  the  original 
Hartwell  algorithm.  Further  test  cases  show  that  the  Moss  algorithm  is  still  fully  eapable  of 
detecting  multiple  eo-channel  emitters,  and,  although  not  conelusive,  indicate  that  it  does 
yield  higher  aceuracy  for  longer  signal  eollection  times. 

Statistieally  eonelusive  results  were  not  achievable  at  this  time  due  to  time  eonstraints 
and  the  high  computational  complexity  involved  in  the  current  implementation  of  the  Moss 
algorithm.  Full  maps,  with  long  signal  lengths,  would  take  approximately  40  days  to  eom- 
pute,  and  even  the  small  and  simple  test  eases  conducted  in  this  work  required  at  least  one 
day  each  to  complete.  This  made  it  unfeasible  to  gather  valid  statistieal  results,  or  to  fully 
investigate  the  long  eollection  times  that  this  algorithm  was  designed  to  use.  So,  although 
all  data  does  indieate  that  improvement  was  achieved,  and  that  the  new  algorithm  merits 
further  investigation,  it  does  require  more  work  before  it  can  be  made  praetieally  useful. 

The  key  limitation  of  the  Moss  algorithm,  as  currently  implemented,  is  the  computation 
time  required  to  find  a  result.  Therefore,  the  first  improvement  neeessary  is  to  make  the 
algorithm  more  effieient,  such  that  the  run-time  is  praetieal.  With  a  reduction  in  run-time 
it  will  be  possible  to  determine  statistical  results  and  to  fully  characterize  the  algorithm. 
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Some  ideas  on  how  to  achieve  this  are: 

1.  Port  the  algorithm  to  a  more  efficient  language  (i.e.  C  or  Fortran). 

2.  Utilize  parallel  processing  [6]. 

3.  Implement  an  adaptive  resolution  regime. 

4.  Investigate  block  updates,  vice  updating  geometry  every  sample. 

Once  this  is  accomplished,  the  reformulated  equations  should  be  re-examined.  The  equa¬ 
tions  formulated  in  this  work  assumed  a  very  basic,  linear  geometry,  and  a  stationary  emit¬ 
ter.  These  assumptions  were  made  for  the  sake  of  mathematical  simplicity  and  clarity 
during  proof  of  concept;  however,  they  are  generally  not  valid  in  real-world  applications. 
Therefore,  the  model  will  need  some  modification  to  more  truly  reflect  reality.  Some  key 
modifications  required  are: 

1.  Expand  the  model  to  three  dimensions. 

2.  Change  the  equations  of  motion  from  linear  to  elliptical  to  account  for  orbital  motion. 

3.  Allow  for  non-zero  emitter  velocities. 

These  modifications  should  not  be  difficult  to  accomplish,  but  will  make  the  mathematical 
model  significantly  more  complex. 

In  conclusion,  it  is  shown  that  the  reformulated  equations  for  TDOA  and  FDOA  are 
accurate,  and  can  be  used  in  a  modified  CAFMAP  algorithm  to  account  for  the  motion  of 
the  collectors  over  time.  The  resulting  Moss  algorithm  is  able  to  geolocate  RF  emitters 
with  accuracy  similar  to  that  of  the  original  Hartwell  algorithm  for  shorter  collection  times, 
and  apparently  with  greater  accuracy  for  long  collection  times.  The  improvement  comes, 
however,  with  the  cost  of  high  computational  complexity.  This  is  significant  enough  that  it 
must  be  mitigated  before  the  new  algorithm  can  be  practically  implemented,  but  the  work 
required  is  merited.  Additionally,  the  approach  used  here,  of  modifying  the  TDOA  and 
FDOA  equations  to  account  for  collector  motion,  should  be  generally  applicable  to  many 
other  methods  of  geolocation  that  estimate  TDOA  and  FDOA  as  a  fundamental  step.  This 
could  lead  to  much  broader  use  of  these  results  once  the  computational  issues  have  been 
solved,  and  improve  the  geolocation  accuracy  of  multiple  methods  in  common  use  today. 
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Figure  4.  CAFMAPS  with  dm  =  100,  SNR  =  10  dB,  and  30  snapshots. 
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CHAPTER  1: 
Introduction 


1.1  Background 

In  today’s  world  of  high-tech  warfare,  we  have  developed  the  ability  to  deploy  virtually 
any  type  of  ordnance  quickly  and  accurately  to  any  location  on  Earth.  We  can  do  this  from 
a  multitude  of  different  platforms,  which  may  be  located  safely  many  miles  away  and  far 
beyond  the  reach  of  retaliation.  These  weapons  have  become  so  accurate  that,  with  caution, 
they  can  be  used  in  the  midst  of  civilian  populations  to  remove  threats  with  minimal  risk 
of  collateral  damage.  What  limits  our  ability  to  neutralize  these  threats  is,  quite  often, 
insufficient  knowledge  of  locations  in  order  to  strike. 

Additionally,  most  conflicts  today  are  asymmetric.  We  are  often  confronting  adversaries 
without  official,  overt  military  forces  and  without  the  large-scale  resources  of  a  nation  state. 
These  covert  groups  substitute  for  their  lack  of  resources  by  moving  swiftly  and  secretly 
behind  the  mask  of  civilian  populations  and  have  repeatedly  proven  themselves  capable  of 
causing  great  harm  to  our  nation  and  world  peace.  While  our  forces  are  more  than  capable 
of  dealing  with  these  people,  they  can  only  do  so  if  they  know  who  and  where  they  are. 

These  facts  indicate  that  one  of  the  most  crucial  elements  of  modem  warfare  is  precisely 
knowing  the  location  of  threats  in  real-time,  or  very  close  to  it.  Many  tools  must  be  used 
in  collaboration  to  achieve  this  knowledge,  one  of  which  is  Signals  Intelligence  (SIGINT). 
SIGINT  methods  involve  collecting  the  Radio  Frequency  (RF)  signals  emitted  by  our  ad¬ 
versary’s  own  use  of  technology  to  gather  information  on  their  location  and  activities.  Of 
particular  importance  to  this  thesis  is  the  method  of  using  satellites  in  orbit  to  collect  emit¬ 
ted  RF  signals  and  to  determine  their  location  on  the  surface  of  the  Earth. 

There  are  many  methods  and  algorithms  for  achieving  this,  most  of  which  use  a  pair  of 
satellites  (collectors)  to  receive  the  signal  and  then  compute  a  Time  Difference  of  Arrival 
(TDOA)  and  Frequency  Difference  of  Arrival  (FDOA)  between  them  [1],  [2],  [3].  These 
parameters  can  then  be  used  to  determine  the  location  of  the  emitter.  This  is  essentially  the 
reverse  of  the  process  used  by  GPS  receivers  on  the  Earth  to  estimate  one’s  position  and 
velocity. 

One  of  these  methods  was  proposed  by  A1  Buczek  [4]  of  the  Naval  Research  Fabora- 
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tory  (NRL)  in  the  early  1990s  and  theoretieally  developed  and  tested  in  a  2005  master’s 
thesis  by  Hartwell  [5]  at  the  Naval  Postgraduate  Sehool  (NPS).  In  this  method,  ealled  the 
Cross  Ambiguity  Function  Geo-Mapping  (CAFMAP)  method  by  Hartwell,  the  Complex 
Ambiguity  Function  (CAF)  of  the  two  signals  received  by  a  pair  of  collectors  is  calculated 
and  then  mapped  to  xy  coordinates  on  the  Earth’s  surface.  The  magnitude  of  the  CAF  forms 
a  mathematical  surface  that  peaks  at  the  correct  TDOA  and  FDOA  values  embedded  in  the 
signals.  When  this  is  mapped  to  the  Earth’s  surface,  the  result  is  a  function  of  position  that 
peaks  at  the  geographical  coordinates  of  the  emitter,  as  demonstrated  in  Eigure  5. 

The  CAFMAP  method  was  successfully  demonstrated  by  Hartwell  [5],  but  it  suffered 
from  a  key  weakness —  its  inability  to  use  long  duration  signals.  Attempts  to  do  so  resulted 
in  a  smearing  of  the  CAF  due  to  a  change  in  the  collector-emitter  geometry  over  the  col¬ 
lection  time.  This  limits  the  best  possible  accuracy  of  Hartwell’s  method  as  the  integration 
time  is  a  fundamental  factor  in  the  precision  of  geolocating  signals  [6] .  Hartwell  was  able 
to  mitigate  this  deficiency,  somewhat,  by  collecting  multiple,  short  snapshots  of  the  RF 
signals  at  spaced-out  intervals,  and  then  non-coherently  summing  the  magnitudes  of  the 
resulting  CAFs.  This  reduced  the  variance  in  the  peaks  of  the  resulting  map  and  provided 
a  better  position  estimate  of  the  emitter,  but  it  was  still  not  an  ideal  solution.  This  thesis 
attempts  to  fundamentally  correct  this  deficiency  by  modifying  the  CAF  algorithm,  such 
that  arbitrarily  long  collection  times  are  made  possible  without  smearing  the  results.  This 
correction  should  yield  a  version  of  the  CAFMAP  method  where  the  estimation  error  can 
be  made  as  small  as  desired  by  simply  increasing  the  collection  duration.  The  trade-off  for 
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Figure  5.  Results  of  Hartwell’s  CAFMAP  algorithm. 
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this  advantage  being  increased  computational  complexity  and  memory  requirements. 


1.2  Objective 

The  primary  objective  of  this  thesis  is  to  derive  a  mathematical  model  of  the  CAP  that 
accounts  for  the  motion  of  the  collectors  and  emitters  in  order  to  prevent  the  smearing  of  the 
CAP  surface  over  long  collection  times.  This  will  allow  the  CAPMAP  method  to  be  used 
with  increased  integration  intervals,  which,  in  turn  will  reduce  the  variance  of  the  TDOA 
and  PDOA  peaks  in  the  resulting  map,  thus  providing  a  more  accurate  geolocation  result. 
This  work  is  focused  on  extending  the  capabilities  of  the  CAPMAP,  but  the  principles 
should  generally  apply  to  any  geolocation  technique  utilizing  the  CAP  to  determine  TDOA 
and  PDOA  values. 

Upon  successfully  deriving  a  mathematical  model,  this  study  implements  it  in  MATLAB 
to  validate  and  characterize  the  results. 

1.3  Related  Work 

Many  papers,  thesis  and  dissertations  have  been  written  on  the  topic  of  CAP  processing  and 
geolocating  RP  signals.  The  most  directly  relevant  work  to  this  thesis  is  that  of  Hartwell, 
where  he  developed  and  demonstrated  the  basic  CAPMAP  algorithm  [5].  This  thesis  orig¬ 
inated  as  an  attempt  to  extend  and  improve  Hartwell’s  work. 

The  fundamental  concepts  and  algorithms  employed  in  geolocation  are  discussed  by 
Loomis  [1],  Price  [2],  and  Groves  [3].  These  fundamentals  have  been  in  practical  use  for 
several  decades  and  still  form  the  foundations  of  SIGINT  geolocation.  The  CAP  is  also 
a  critical  concept  in  this  thesis  and  its  definition,  uses  and  methods  of  computation  are 
developed  in  [6],  [7],  [8],  and  [9].  Skolnik  also  provides  some  useful  insight  into  the  uses 
of  the  CAP  in  RADAR  [10],  where  the  function  found  its  original  use. 

Johnson’s  Naval  Postgraduate  School  thesis  [11]  is  also  of  great  value.  Johnson  cre¬ 
ated  and  tested  the  signal  generating  functions  used  in  this  work  and  provided  a  practical 
implementation  of  using  the  CAP  to  perform  geolocation. 
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1.4  Thesis  Organization 

This  thesis  is  organized  into  five  chapters.  Chapter  2  gives  an  overview  of  the  CAP  and 
a  demonstration  of  how  it  can  be  used  to  estimate  TDOA  and  FDOA  parameters  simul¬ 
taneously.  It  also  introduces  the  basic  factors  affecting  the  accuracy  of  the  CAP.  Chap¬ 
ter  3  develops  the  mathematical  model  for  the  new  CAFMAP  algorithm  and  discusses  its 
implementation  in  MATLAB.  Chapter  4  presents  results  and  compares  them  to  those  of 
Hartwell’s  original  method.  Chapter  5  concludes  this  thesis  by  presenting  some  ideas  for 
future  work  and  ways  that  could  make  this  new  algorithm  more  useful  in  a  practical  appli¬ 
cation. 
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CHAPTER  2: 

Overview  of  the  CAP 

This  study  is  focused  on  extending  the  capabilities  of  the  CAP.  Thus,  the  first  step  is 
to  arrive  at  a  good  understanding  of  what  a  CAP  is  and  how  it  is  used  in  the  field  of 
geolocation.  In  essence  the  CAP  is  a  correlation  function  that  takes  a  time  and  frequency 
offset  as  its  independent  variables  and  peaks  when  the  values  of  those  variables  match  the 
inherent  offsets  embedded  in  a  received  signal.  It  has  found  widespread  use  in  the  fields 
of  radar  detection  and  geolocation  because  of  its  ability  to  simultaneously  estimate  both 
TDOA  and  PDOA  values  [1],  [6],  [7],  [9],  [10],  [11]. 

2.1  Definition 

There  are  multiple  definitions  of  the  CAP  but  the  most  common,  and  the  one  that  will  be 
used  in  this  text,  is  that  given  by  Stein  [6]  as 

A(t,/)=  r  +  (2.1) 

Jo 

The  integration  limit,  T,  is  the  length  of  time  over  which  the  signals  are  collected.  Purther- 
more,  when  applied  specifically  to  geolocating  a  received  signal: 

51  (t)  =  analytic  signal  received  at  collector  1, 

S2{t)  =  analytic  signal  received  at  collector  2, 

T  =  time  difference,  and 
/  =  frequency  offset. 
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2.2  Simple  Example 

To  understand  how  this  works  we  will  use  a  simple  signal,  reet(t),  as  an  example.  This 
signal  is  shown  in  Figure  6  and  is  defined  as: 


reet(t) 


1  if— -<t<- 

Ai  2—^  —  2 

0,  Otherwise. 


(2.2) 


In  a  geoloeation  problem,  a  signal  is  emitted  by  a  transmitter.  This  signal  is  then  deteeted  by 
two  separate  eolleetors  with  different  positions  and/or  veloeities.  These  reeeived  signals  are 
5i(t)  and  S2{t).  Both  reeeived  signals  are  the  original  reot(t)  signal  that  was  transmitted, 
but  they  will  have  slightly  different  times  of  arrival  and  frequeneies  in  addition  to  some 
added  noise.  The  times  of  arrival  are  different  beeause  the  two  eolleetors  are  in  different 
positions,  therefore  the  path  length  from  the  emitter  to  eaeh  eolleetor  is  different.  The 
differenee  in  path  length  eauses  a  differenee  in  time  of  arrival.  Similarly,  the  eolleetors  will 
typieally  have  different  veloeity  veetors.  The  relative  veloeity  between  eaeh  eolleetor  and 
the  emitter  will  eause  a  Doppler  shift,  but  sinee  the  relative  veloeity  veetors  are  not  equal, 
the  amount  of  Doppler  shift  will  not  be  either.  This  leads  to  a  differenee  in  frequeney 
between  the  two  reeeived  signals.  The  following  example  shows  how  these  offsets  are 
determined  using  the  CAF.  For  this  example,  the  following  terms  will  be  used: 


T  =  possible  time  differenee  parameters, 

/  =  possible  frequeney  offset  parameters. 

At  =  aetual  time  differenee  between  signals,  and 
A/  =  aetual  frequeney  differenee  between  signals. 


The  objeetive  is  to  ealeulate  the  CAF,  whieh  is  defined  in  (2.1).  This  definition  follows 
that  given  by  Stein  [6],  and  is  used  in  this  work,  but  the  funetion  ean  also  be  deseribed  in 
a  more  general  form,  with  infinite  bounds  of  integration,  sueh  as  the  definition  found  in 
Ramirez  et  al.  [9].  The  two  definitions  are  equivalent,  in  praetiee,  as  the  eolleeted  signals 
are  of  finite  length,  starting  at  t  =  0  when  the  eolleetion  starts,  and  ending  att  =  T  when  the 
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rect(t) 


_ 1 _ 

< - > 

-1/2  i  1/2  ^ 

Figure  6.  rect(f)  Signal. 

collection  ends.  In  practical  applications,  the  Stein  definition  makes  more  intuitive  sense, 
and  is  the  definition  we  use;  but  in  this  example,  we  will  revert  to  the  infinite  bounds  for 
mathematical  convenience.  Understand  that  they  are  equivalent  in  practice.  Thus,  we  let 


/CX3 

+  (2.3) 

-oo 

where 

5i(t)  =  and  (2.4) 

Slit)  =  rect(t  —  (2.5) 

For  simplicity,  the  problem  will  be  broken  down  into  segments.  Let, 

mit)  =  Slit)  slit +  t),  (2.6) 

which  is  known  as  the  mixing  function,  and  yields 

mit)  =  siit)slit  +  z)  (2.7) 

=  rect(t)e'^^^'^‘'^rect(t  -At  +  (2.8) 

=  rect(t)rect(t  -  At  -f  (2.9) 


The  mixing  function  is  then  combined  with  the  exponential  term  of  the  CAF,  taking  the 
frequency  offset  into  account,  yielding 
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=  rect(Orect(?  -  A?  +  (2.10) 

=  rect(?)rect(?  -  A?  +  (2.1 1) 

The  integral  is  then  calculated  yielding 

/CX3 

rect(t)rect(t  —  At  +  (2.12) 

-oo 

/oo 

rect(t)rect(t  —  (At  —  (2.13) 

-oo 

In  order  to  develop  an  understanding  of  the  CAP  it  is  useful  to  look  at  the  parts  of 
this  integral,  rather  than  just  finding  the  analytic  solution.  To  see  the  results  of  the  above 
integration  we  first  neglect  the  frequency  shift  and  just  examine  the  Fourier  transform  of 
the  two  rect()  functions. 

If  |At  —  t|  >  1  then, 


rect(t)rect(t  —  (At  —  t))  =0  (2.14) 

^A(t,/)=0.  (2.15) 

This  is  illustrated  in  Figure  7.  If  |At  —  t|  >  1  then  the  two  rect()  functions  have  zero 
overlap  and  their  product  is  zero,  making  the  CAP  equal  zero. 


rect(t) 

^  rect(t-l) 

-1 

n 

1/ 

f2  ^ 

Figure  7.  rect(t)rect(t  —  1). 
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Otherwise,  if  |At  —  t|  =0  then, 

reet(t)reet  ((t  —  (At  —  t))  =  reet^(t) 

=  reet(t). 


and  sinee 


noo 

/  XQCi{t)e^  dt  =  sine(/  —  A/) 

J  oo 

A{z  =  At,/)  =  — A/). 


(2.16) 

(2.17) 


(2.18) 

(2.19) 


The  first  faetor  is  a  eomplex  phasor,  meaning  that  it  will  ehange  the  phase  of  the  result,  but 
the  magnitude  will  remain  unity.  The  seeond  faetor  is  the  well-known  sinc()  funetion,  and 
its  magnitude  will  equal  one  if  \  f  —  Af\  =  0,  and  will  oseillate  with  deereasing  magnitude 
otherwise.  From  this,  it  is  seen  that  |A(t,/)|  will  always  be  maximized  when  both  T  =  At 
and/  =  A/. 

The  final  ease  is  where  0  <  |At  —  t|  <  1.  As  shown  in  Figure  8,  the  mixing  funetion 
m(t)  =  reot(t)reet  (t  —  (At  —  t))  beeomes  the  shifted  and  sealed  reet(t)  function 


m(t)  =  rect 


At-T 


1  -  |At-T| 


(2.20) 


This  gives  a  CAF  equation  of 

pcxy  _ 

A  (0  <  |At  —  t|  <  1,/)  =  — 

J  oo  \  i 


Af-T 


|At-T| 


^j2Tt{f-Af)t^^_  (2.21) 


It  can  be  seen  that  the  integral  in  (2.21)  is  just  the  Fourier  Transform  (FT)  of  (2.20),  from 
which  we  can  make  use  of  the  FT  shift  and  scaling  properties.  Doing  so  yields 


rect 


t- 


Ar-T 


1  -  |At-T| 


e 


A/)|^  _  |At-T||sinc{[l  -  |At-T|]  (/-A/)} 

(2.22) 
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Figures.  rect(f)rect(f  —  (Af  —  r)). 


^  A(0  <  |A?-  t|  <  1,/)  = 

^-j2K{f,-Af)T^-j2K{^){f-Af)  1 1  _  1^^  _  ^1  |sinc  { [1  -  I  A?  -  t|]  (/  -  A/) } .  (2.23) 

Importantly,  the  boundaries  of  (2.23)  agree  with  the  previous  oases.  When  |At  —  t|  =  0, 
then  the  term  [1  — |Ar  — t|]  equals  1,  and  the  sine  term  becomes  sine  (/  —  A/).  This  is 
maximized  when  \  f  —  Af\  =  0.  When  |At  —  t|  =  1,  then  [1  —  |At  —  t|]  =  0  and  forces  the 
entire  function  to  zero.  In  between  these  boundaries,  it  can  be  seen  that  as  |  At  —  t|  increases, 
the  coefficient  term  [1  —  |Ar  —  t|]  becomes  smaller,  which  reduces  the  magnitude  of  the 
CAP.  Therefore,  the  CAP  will  have  a  higher  magnitude  when  the  difference  between  T 
and  At  is  smallest.  At  the  same  time,  the  sinc()  function  is  maximized  when  the  difference 
between  /  and  A/  is  smallest  and  decreases  as  the  difference  is  increased.  Additionally, 
both  exponential  factors  will  cause  oscillations,  but  their  magnitude  is  1  when  |Af  —  t|  =  0 
and  1/  —  A/I  =0  and  can  never  be  greater  for  any  other  values. 

So  it  can  be  seen  in  this  example  that  when  the  CAP  surface  is  plotted  with  respect  to  T 
and  /,  as  shown  in  Pigure  9,  there  will  be  a  peak  at  the  location  corresponding  to  the  actual 
signal  offsets.  At  and  A/.  This  demonstrates  that  by  numerically  calculating  the  CAP  of  a 
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signal,  received  by  two  different  collectors,  we  can  arrive  at  a  good  estimate  of  the  values 
for  both  the  TDOA  and  the  FDOA.  There  will,  however,  be  noise  is  this  estimate,  mostly 
due  to  noise  inherent  in  the  signals  themselves.  Therefore,  one  of  our  goals  when  utilizing 
the  CAF  for  geolocation  is  to  minimize  the  variance  in  the  TDOA  and  FDOA  estimates, 
in  order  to  maximize  accuracy.  Of  note  in  Figure  9  is  that  the  variance  in  the  FDOA  is 
practically  zero.  This  is  due  to  the  infinite  bandwidth  of  the  rect(0  signals  used  in  this 
example  and  this  relationship  will  be  further  discussed  in  the  next  section. 

2.3  CAF  Performance 

Along  with  his  work  on  how  to  compute  CAF  surfaces,  Stein  also  derived  the  theoretical 
accuracies  of  both  TDOA  and  FDOA  estimates.  These  are 


<yTDOA 

(2.24) 

Tg  a/BTy' 

(2.25) 

with  the  terms  defined  as: 


T  =  signal  collection  time,  i.e.  length  of  5i(t)  and  52(t), 

B  =  noise  bandwidth  at  receiver  input,  assumed  the  same  for  both  collectors, 
j8  =  RMS  bandwidth  of  the  received  signal, 

Tg  =  RMS  integration  time,  and 
7  =  effective  input  SNR. 

Furthermore, 


1 

2 

(2.26) 


where  Ws(/)  =  signal  power  density  spectrum. 
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X:  0.104 
.  Y:  -10.01 
Z:  1237 


FDOA  (Hertz)  TDOA  (Seconds) 


Figure  9.  CAF  of  the  example  rect(f)  signals,  with  offsets  of  10  Hz  in  frequency,  and  0.1  seconds 
in  time. 


T.  =  In 


roo\uit)\^dt 


where  u{t)  is  the  waveform  that  was  originally  transmitted,  and 


(2.27) 


1 

7 


1  r  1 


1 

- h 

72 


1 

7172J  ’ 


(2.28) 


where  71  and  72  are  the  signal  to  noise  ratio’s  in  the  signal  reeeived  by  eaeh  eolleetor  [6]. 

It  is  clear  that  in  order  to  minimize  the  variance  in  both  TDOA  and  FDOA  dimensions, 
we  need  to  maximize  the  five  terms  in  (2.24)  and  (2.25).  B,  /3  and  7  are  all  dependent  on  the 
signal  that  is  transmitted  and,  therefore,  out  of  the  collectors’  control.  That  leaves  T  and 
Tg,  which  are  both  measures  of  the  signal  collection  time.  Thus,  if  we  want  to  reduce  the 
variance  in  the  CAF  peaks,  and  thereby  improve  the  geolocation  estimate,  it  is  necessary  to 
increase  the  signal  collection  time  as  much  as  possible. 

This  chapter  has  developed  an  understanding  of  the  CAF  and  the  key  factors  affecting  its 
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accuracy.  The  next  chapter  will  introduce  Hartwell’s  CAFMAP  method  of  geolocation,  and 
explain  its  inherent  limitations  regarding  signal  collection  time.  It  will  then  explore  how 
we  might  reformulate  the  CAF,  such  that  the  inherent  motion  of  the  collectors  is  accounted 
for.  In  principle,  this  will  allow  us  to  collect  signals  for  as  long  as  necessary  to  achieve 
our  desired  level  of  accuracy.  In  reality  there  will  be  other  limiting  factors,  such  as  latency, 
memory,  and  processing  capability,  but  the  physical  limitations  of  the  CAFMAP  method 
will  have  been  overcome. 
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CHAPTER  3: 

A  New  CAFMAP  Algorithm 


Hartwell  demonstrated  a  novel  method  of  using  the  CAP  for  geoloeation  in  his  master’s 
thesis  whieh  he  ealled  the  CAFMAP  [5].  Traditional  methods  use  the  CAP  to  determine  the 
TDOA  and  FDOA  values,  and  then  utilize  a  follow-on  method,  sueh  as  Newton-Raphson,  to 
invert  the  TDOA  and  FDOA  equations  in  order  to  arrive  at  geographic  coordinates  [1],  [2]. 
The  CAFMAP  method,  instead,  maps  the  CAP  directly  to  geographic  coordinates.  This 
results  in  a  surface,  or  image,  that  peaks  at  the  most  likely  location  of  the  RF  emitter. 

Hartwell’s  approach  for  doing  this  was  to  create  a  grid  map  of  x  and  y  geographic  coor¬ 
dinates  and  calculate  what  the  expected  TDOA  and  FDOA  values  would  be  if  the  emitter 
were  located  in  each  of  those  grid  positions.  He  then  calculated  the  actual  CAP  surface 
from  the  collected  signals  and,  by  cross-referencing  the  TDOA  and  FDOA  values  of  this 
surface  and  the  grid  points,  was  able  to  map  the  CAP  into  x  and  y  coordinates.  The  result 
was  the  CAP  surface  on  an  xy-coordinate,  geographic  plane.  Hartwell’s  results  were  quite 
promising,  particularly  when  dealing  with  multiple  co-channel  emitters  as  shown  in  Figure 
10. 

Unfortunately,  it  also  relied  upon  constant  signal  parameters  to  be  able  to  perform  the 
cross-referencing  and  also  to  be  able  to  calculate  the  CAP  efficiently  using  Fast  Fourier 
Transform  (FFT)  techniques.  The  implications  of  this  restriction  are  shown  in  Figure  11. 
As  the  signal  lengths  are  increased  the  CAP  peaks  become  greater  in  magnitude  and  drop 
off  more  sharply.  This  means  that  the  geoloeation  accuracy  will  be  improved,  as  expected. 


'  4 
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Figure  10.  Results  of  Hartwell’s  CAFMAP  algorithm. 
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for  longer  collection  and  integration  times.  However,  when  the  signal  lengths  reach  ap¬ 
proximately  2^®  samples  long  we  can  begin  to  see  more  noise  in  the  form  of  spurious 
peaks  away  from  the  true  emitter  location.  A  signal  length  of  2^^  samples  corresponds  to 
a  collection  time  of  approximately  10.5  seconds  for  the  typical  sampling  rate  of  100  kHz 
used  in  this  example.  This  noise  arises  due  to  the  fact  that  the  collector-emitter  geometry 
has  changed  over  the  length  of  the  signal  collection  interval,  and  it  results  in  geolocation 
accuracy  becoming  worse  with  longer  collection  times,  rather  than  better. 

Hartwell  recognized  this  limitation  and  made  an  effort  to  mitigate  it.  He  did  this  by 
breaking  the  signal  collection  up  into  short  duration  snapshots,  with  a  time  gap  between 
each.  He  then  calculated  the  CAFMAP  for  each  of  these  snapshots  and  non-coherently 
summed  their  magnitudes.  The  idea  of  this  being  that  each  snapshot  would  have  its  own 
noise,  but  the  same  peaks  where  emitters  were  located.  Thus,  by  summing  the  snapshot 
magnitudes,  the  surface  energy  at  the  true  emitter  locations  would  keep  increasing  at  a 
greater  rate  than  the  side  lobes.  This  mitigation  technique  was  partially  successful.  It  al¬ 
lowed  for  longer  effective  signal  collections,  and  yielded  correspondingly  narrower  peaks, 
but  it  also  tended  to  create  its  own  noise.  Much  of  this  is  due  to  the  non-coherent  summing 
used  by  Hartwell.  With  no  phase  information  there  is  no  destructive  interference,  so  the 
noise  floor  is  always  increasing.  Eventually  this  makes  it  difficult  to  distinguish  peaks  from 
the  noise.  Additionally  there  is  a  limit  to  how  narrow  the  peaks  can  be  made  with  this 
technique.  Since  each  CAP  snapshot  must  be  from  a  short  duration  signal,  there  will  be  a 
relatively  greater  variance  in  the  peaks  being  summed.  Summing  them  will  cause  the  center 
of  the  peak  to  rise  faster  than  the  surrounding  surface,  but  it  will  not  remove  the  energy  of 
the  individual  wide  peaks  [5]. 

The  goal  of  this  thesis  is  to  completely  remove  this  physical  limitation.  This  is  accom¬ 
plished  by  modifying  the  CAP  equations  used  in  the  CAPMAP  method,  such  that  system 
geometry  is  updated  throughout  the  signal  collection  period.  In  principle,  this  allows  us  to 
collect  signals  for  as  long  as  necessary  to  achieve  our  desired  level  of  accuracy.  In  reality 
there  will  be  other  limiting  factors,  such  as  latency,  memory,  and  processing  capability,  but 
the  physical  limitations  of  the  CAPMAP  method  will  have  been  overcome. 

Initially,  a  variety  of  modifications  to  Hartwell’s  established  algorithm  were  attempted 
to  achieve  the  desired  result.  After  some  effort,  however,  it  was  realized  that  the  basic 
architecture  of  Hartwell’s  algorithm  was  incompatible  with  the  objective.  Therefore,  the 
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basic  formulation  of  the  problem  was  revisited  as  follows. 


3.1  Developing  the  Basic  Model 

The  critical  aspect  of  the  problem  is  the  geometry  between  the  collectors  and  the  emitter.  A 
generic  example  geometry  is  shown  in  Figure  12,  which  will  be  used  to  derive  a  new  model. 
A  three-dimensional  (3D)  geometry  is  used  for  generality,  however,  certain  restrictions 
will  be  enacted  during  this  development  in  order  to  maintain  mathematical  simplicity.  The 
standard  CAF  is  a  function  of  TDOA  and  FDOA,  but  our  objective  requires  that  it  be  made 
into  a  function  of  position  and  time,  or,  x,  y,  z  and  t.  Accomplishing  this  is  the  first  step  in 
developing  the  new  model. 
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Figure  11.  Hartwell  CAFMAP  results  for  varying  signal  lengths. 
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3.1.1  TDOA 

For  simplicity,  TDOA  and  FDOA  will  be  transformed  separately,  starting  with  TDOA. 
TDOA  is  simply  the  difference  in  time  that  is  required  for  the  signal  to  travel  from  the 
emitter  to  the  two  collectors.  That  is. 


T  = 


K2  -  n 


(3.1) 


where  c  =  speed  of  light  2.998  x  IQ^m/s  [1],  [2]. 

Or,  in  the  terms  of  Figure  12  and  using  the  velocity  of  the  collectors  and  emitter,  which 
are  allowed  to  be  time- variant. 


T(t) 


Pe{t)-Pc,{t)  -  J  P,[t)  -  Pc,{t) 


Pe-\-Ve't]  —  \  Pc2  +^02'^ 


(3.2) 

(3.3) 


For  the  remainder  of  this  work  it  is  assumed  that  Vg  =  0,  i.e.  only  stationary  emitters  are 
considered.  Examining  the  effects  of  this  assumption  will  be  left  as  a  future  work.  With 
this  assumption  (3.2)  simplifies  back  to 


where 


^  1^2(01 -In (01 


n  (0  =Pe-  (Pci  +  Vci  ■  and 
n(0  =Pe-  (Pc2+V^ci  . 


(3.4) 


(3.5) 

(3.6) 


This  successfully  converts  the  TDOA  into  a  function  of  position  and  time.  The  velocities 
of  the  collectors  must  be  known,  but  in  general  they  are  and  this  is  not  a  problem. 
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Figure  12.  Example  geolocation  geometry. 


3.1.2  FDOA 

As  stated  by  Loomis  [1],  the  Doppler  shift  in  signal  frequency,  A^,  due  to  the  motion  of  a 
collector  is  given  by 


A(^  = 


fcV-r 


(3.7) 


where  fc  is  the  carrier  frequency  of  the  emitted  signal,  and  c  is  the  speed  of  light.  In  the 
terms  of  Figure  12,  this  means  that  the  frequency  shifts  for  each  collected  signal  are 


A  A  fciVci-ri 

A01  =  1  ,  and 

(3.8) 

c  Ini 

(3.9) 

c  |r2| 
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Since  the  two  reeeived  signals  are  simply  different  eopies  of  the  same  original  signal,  fc^  — 
fc2  —  fc-  Thus,  the  frequeneies  of  the  signals  reeeived  by  the  two  eolleetors  are 

^i=fc  +  Ah  =fc  +  -  and  (3.10) 

c  \ri  \ 

(l)2  =  fc  +  A(l>2=fc+-'^^.  (3.11) 

c  \r2\ 

The  FDOA  is  simply  the  difference  in  these  two  reeeived  signal  frequeneies,  and  therefore 


f  =  h-h  =  - 

c 


r  ■  ■ 

ri  Vci  ■  n  ] 

[ 

1^2! 

In 

j 

(3.12) 


When  n  and  r2  are  allowed  to  be  time-variant,  and  are  defined  by  (3.5)  and  (3.6),  (3.12) 
becomes 

r ■  rn(t)  v2.  ■  1 

(3.13) 


fjL 

'vh-nit) 

vci  ■  h  (t) ' 

c 

.  ln(0l 

In  (01  . 

Thus  we  have  aeeomplished  our  goal,  as  (3.13)  now  provides  FDOA  in  terms  of  x,  y,  z  and 
t. 


3.1.3  Modified  CAF 

With  (3.2)  and  (3.13)  for  TDOA  and  FDOA,  respeetively,  it  is  now  possible  to  determine  a 
form  of  the  CAF  with  position  as  its  independent  variables.  Starting  from  the  general  form 
of  the  CAF  as  given  in  (2.1), 

A(t,/)=  +  (3.14) 

Jo 

where  the  signals  (t)  and  52 (t)  are  of  length  T,  thus  bounding  the  integral. 

To  be  praetieally  implemented  (3.14)  must  be  eonverted  to  a  diserete  form  to  allow  for 
digital  signal  analysis  and  eomputation.  With  fs  being  the  sampling  frequeney  of  the  eol- 
leeted  signals,  and  Ts  =  V the  sample  duration,  let 


t  =  nTs 


(3.15) 
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and 


(3.16) 


f=^A 
J  N  ’ 

where  N  is  the  number  of  collected  samples  for  each  signal,  such  that  T  =  NTs,  and  k  is  the 
digital  FDOA.  These  substitutions  lead  to  the  discrete  CAF, 


N~l 

A{x,k)  =  ^  5i  [n]s^ 

n=0 


X 

n-\ - 

Ts 


(3.17) 


Substituting  in  (3.2)  and  (3.13),  letting  hq  equal  the  first  sample  of  the  collected  signals, 
and  noting  that  N  =  fsT,  yields  the  modified  CAF  function. 


A{x,y)  = 

Vc2-^(nTs)  Vej-r*[(n7i) 
l''2(«rj)|  kl(n7i)| 

(3.18) 

In  (3.18)  the  term  j  must  be  rounded  to  the  nearest  integer  because  it  is  in 

the  digital  domain  and  corresponds  to  an  offset  in  the  sample  index,  meaning  that  fractional 
offsets  are  not  allowed. 

3.2  Mapping  Implementation 

Equation  (3.18)  is  implemented  in  two  functions,  LongIntCafMap.m  and  caf_map.m, 
which  can  be  found  in  the  Appendix.  These  functions  will  collectively  be  referred  to  as 
the  Moss  algorithm  for  the  remainder  of  this  work.  The  LongIntCafMap.m  function  de¬ 
fines  the  geometry  of  the  collectors  and  emitter,  the  parameters  of  the  signal,  and  the  length 
and  number  of  snapshots  that  will  be  collected.  It  then  generates  the  "received"  signals 
for  each  collector  using  the  sig_gen.m  function,  written  by  Johnson  in  his  thesis  on  CAF 
processing  [11].  Although,  in  principle,  the  Moss  algorithm  is  able  to  collect  one  long  sig¬ 
nal,  this  implementation  continues  to  break  it  up  into  smaller  snapshots.  This  is  done  for 
two  reasons.  First,  we  want  to  compare  the  Hartwell  and  Moss  algorithms  to  determine  any 
improvement  made.  Thus,  we  want  to  keep  as  many  parameters  as  possible  the  same.  If  the 
Moss  algorithm  shows  improvement  using  snapshots,  then  we  can  be  reasonably  confident 


no+fsT-l 

£  si[n]s*2 


n-hROUND  ki(n7;)|\ 

e 

V  cTs  _ 
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that  it  would  show  even  greater  improvement  with  just  one,  long  eolleetion.  Seeond,  the 
trade-off  for  updating  the  system  geometry  every  sample  is  very  high  eomputational  eom- 
plexity.  This  means  that  one  long  eolleetion  would  take  prohibitively  long  to  simulate,  and 
eontinuing  to  use  the  snapshots  cuts  down  on  this  processing  time  greatly.  After  generating 
the  "received"  snapshots,  the  LongIntCafMap.m  function  passes  them  to  the  caf_map.m 
function. 

The  caf_map.m  function  accepts  the  collected  signals  and  parameters  for  the  size  and 
resolution  of  the  map.  It  then  loops  through  every  x  and  y  position  in  the  grid  and  calcu¬ 
lates  the  CAP  value  at  that  position  using  (3.18).  Unfortunately,  this  must  be  done  in  a 
loop  rather  than  utilizing  the  FFT  or  cross-correlation  methods  because  the  parameters  are 
changing  in  time,  meaning  that  they  must  be  updated  every  sample.  As  a  result,  this  algo¬ 
rithm  is  unable  to  take  advantage  of  the  computational  efficiency  of  those  algorithms  and 
must  utilize  brute  force,  as  it  were,  by  looping  through  every  signal  sample  and  computing 
each  element  individually.  This  causes  a  very  high  level  of  computational  complexity,  and 
thus,  the  algorithm  is  very  slow.  This  is  further  discussed  in  the  succeeding  chapters.  After 
calculating  the  CAP  for  each  snapshot,  caf_map.m  returns  it  to  LongIntCafMap.m,  which 
coherently  sums  the  values  of  all  snapshots  to  yield  the  final  CAFMAP  result. 

3.3  Key  Differences 

The  next  chapter  will  discuss  the  Moss  algorithm  in  more  detail  and  compare  its  results 
to  those  of  the  Hartwell  algorithm.  However,  it  is  important  to  note  the  key  differences 
between  the  two  before  proceeding. 

1.  Hartwell  algorithm: 

(a)  Snapshots  are  required. 

(b)  Snapshots  are  summed  non-coherently. 

(c)  The  CAP  surface  is  binned  into  defined  grid  locations. 

(d)  The  FFT  is  used  for  computing  CAFs;  very  efficient. 

2.  Moss  algorithm: 

(a)  Snapshots  are  not  required;  single,  long  duration  samples  may  be  used  instead. 

(b)  When  snapshots  are  used,  they  are  coherently  summed  to  maintain  phase  infor¬ 
mation. 

(c)  The  CAP  surface  is  computed  directly  in  the  xy  domain. 
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(d)  Each  element  of  the  CAP  must  be  computed  individually  in  loops;  high  com¬ 
putational  cost. 

In  short,  the  Moss  algorithm  gains  accuracy  and  resolution  at  the  cost  of  computational 
complexity.  The  implications  of  this  trade-off  are  fully  examined  in  the  next  chapter. 
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CHAPTER  4: 
Results 


The  initial  results  show  that  the  Moss  algorithm  does  indeed  locate  RF  emitters  and  main¬ 
tains  the  CAFMAP’s  ability  to  simultaneously  detect  co-channel  emitters  as  shown  in  Fig¬ 
ure  13.  This  justifies  the  direct  mapping  solution  used  and  opens  the  door  for  further  work 
on  the  algorithm. 

Having  verified  that  the  algorithm  works  as  intended,  the  first  task  is  to  characterize  it 
and  check  its  performance  compared  to  the  original  Hartwell  algorithm.  The  computational 
complexity  of  the  Moss  algorithm  turns  out  to  have  a  large  bearing  on  the  ability  to  char¬ 
acterize  the  accuracy,  so  it  is  discussed  first.  This  is  followed  by  a  side-by-side  comparison 
of  the  results,  obtained  by  using  the  two  algorithms  and  a  discussion  on  their  differences. 

4.1  Computational  Complexity 

Listing  4.1  shows  the  heart  of  the  Moss  algorithm  as  currently  implemented  in  MATLAB. 
In  order  to  analyze  the  computational  complexity  of  the  algorithm  the  methodology  chosen 
is  as  follows.  We  assume  that  each  operation  requires  one  unit  of  time  to  perform,  mean¬ 
ing  that  addition,  subtraction,  multiplication,  division,  assignments,  etc.  all  take  the  same 
amount  of  time  to  perform.  How  true  this  assumption  is  will  vary  depending  on  the  proces¬ 
sor  being  used,  but  for  an  initial  approximation,  it  is  sufficient  and  simplifies  the  analysis. 
Additionally,  the  more  complex  operations,  such  as  norms  and  vector  multiplies,  are  bro¬ 
ken  up  into  their  component  basic  operations.  For  example,  the  norm  of  a  three  element 
vector  is  norm(x)  =  \Jx\+x\  +x'^,  so  it  is  counted  as  three  multiply,  two  add,  and  one  root 
operations. 

This  done,  the  number  of  operations  performed  in  each  loop  is  counted  and  the  result 
multiplied  by  the  number  of  times  that  the  loop  is  executed.  The  results  of  this  are  shown 
in  Table  1.  When  the  number  of  all  operations  are  summed  we  arrive  at  #  Ops  =  1  4- 
NxNy{4  +  78A),  where  Nx  and  Ny  are  the  number  of  gridpoints  in  the  xy  grid,  and  N  is  the 
number  of  collected  samples  for  the  signals.  As  an  example,  if  we  assume  that  the  area 
of  interest  is  an  area  50  km  wide  and  20  km  deep  at  a  resolution  of  10  meters,  and  the 
collection  time  is  approximately  40  seconds  at  100  kHz  sampling  rate,  then  the  parameters 
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would  be: 


=5000, 

Ny  =2000  and 
N  =2^2  _ 

The  total  number  of  operations  to  be  performed  is  1  +  NxNy{A  +  7SN)  =  3.27  x  10^^  and  if 
it  is  assumed  that  one  operation  requires  one  nanosecond  to  perform  then  the  total  run  time 
will  be  (3.27  x  10^^)  (1  x  10^^)  =  3.27  x  10^  seconds,  or  37.86  days.  This  is  obviously 
impractical,  and  time  was  not  available  to  experimentally  confirm  this  run-time.  However, 
experiments  were  performed  with  either  lower  resolutions  or  shorter  sample  lengths  and 
their  run-times  in  MATLAB  are  shown  in  Figures  14  and  15.  These  figures  show  the 
run-time,  in  seconds,  of  each  function  called  during  the  simulation.  LonglntCafMap  is  the 
main  function,  and  its  total  time  is  the  total  time  taken  for  the  entire  simulation.  The  results, 
shown  in  these  figures,  confirm  the  complexity  of  the  Moss  algorithm  and  demonstrate  its 
impracticality  as  currently  implemented. 
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Listing  4.1.  Moss  algorithm  in  MATLAB  code. 


for  ii  =  1 : Nx 

for  j  j  =  1 : Ny 

m  =  zeros { 1 , N) ; 
exponent  =  zeros (1,N); 

Pe  =  [ indexX ( ii ) , indexY ( j j )  ,  0 ]  ; 

S2_conj  =  conj(S2); 
for  nn  =  1:  length (S2) 

rl  =  Pe  —  (Pcl+Vcl . *nn*Ts ) ; 
r2  =  Pe  —  (Pc2+Vc2 . *nn*Ts )  ; 

TauValue(nn)  =  round ( (norm { r2) —norm (rl )) /c) ; 
if  TauValue (nn)  >=  0  &&  (nn+TauValue (nn) )  <  N 

m(nn)  =  SI  (nn)  *  S2_con j (nn+TauValue (nn) )  ; 
elseif  (nn— TauValue (nn)  )  <  N 

m(nn)  =  SI (nn— TauValue (nn) )  *  S2_conj (nn) ; 

end 

FDOA(nn)  =  Fo/c* (Vc2 *r2 ' /norm (r2) —Vcl*rl ' /norm (rl ))  ; 
exponent (nn)  =  exp (— 1 j*2*pi*FD0A (nn) /Fs* (nn— 1) )  ; 

end 

product  =  m. *exponent; 

CAF(ii,jj)  =  sum (product ) ; 

end 

end 


Table  1.  Number  of  operations  performed. 


Operation 

#  of  times  performed 

Assign 

NxNy{\0N  +  3) 

Add/Sub 

NxNy29N 

Multiply 

NxNy{26N  +  N+l) 

Divide 

NxNySN 

Sqrt 

NxNyAN 

Compare 

NxNy2N 

Exponent 

NxNyN 

Nx  =  #ofX  grid  positions  mapped. 
=  #  of  Y  grid  positions  mapped. 
N  =  #  of  signal  samples  taken. 
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;  Start  Profiling  Run  this  code;  ■»  ®  Profile  time  75850  sec 


Profile  Summary 

Generated  20-May-2014  11:00:16  using  cpu  time. 


Function  Name 

Calls 

Total  Time 

Self  Time* 

Total  Time  Plot 
(dark  band  =  self  time) 

LonaintCaflvIap 

1 

75847.664  s 

0.081  s 

caf  map 

10 

75843.880  s 

75843.880  s 

sia  aen 

10 

3.572  s 

3.543  s 

mesh 

1 

0.052  s 

0.002  s 

close 

1 

0.049  s 

0.001  s 

setdiff 

4 

0.040  s 

0.006  s 

sefdiff>setdiffleaacv 

4 

0.034  s 

0.015  s 

newolot 

1 

0.031  s 

0.002  s 

Figure  14.  N  =  2^^,  dm  =  100  m,  full  map. 


^  Profiler 


iHl 
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Figure  15.  N  =  2^^,  dm  =  100  m,  small  map. 


28 


4.2  Moss  Algorithm  versus  Hartwell  Algorithm 

As  discussed  in  the  previous  section,  a  full  characterization  of  the  Moss  algorithm  was 
unfeasible.  In  lieu  of  this,  a  selection  of  test  cases  were  run  with  both  the  Moss  and  the 
Hartwell  algorithms  so  that  a  rough  comparison  might  be  made  and  the  general  charac¬ 
teristics  of  the  Moss  algorithm  determined.  These  tests  were  conducted  with  a  lead-trail 
collector  geometry,  with  the  geometry  parameters  defined  in  Table  2.  The  transmitted  sig¬ 
nal  is  Binary  Phase  Shift  Keying  (BPSK)  modulated  with  parameters  as  shown  in  Table  3. 
All  tests  are  performed  using  both  the  re-created  Hartwell  algorithm  and  the  Moss  algo¬ 
rithm,  with  identical  parameters,  such  that  a  direct  comparison  can  be  made. 


Table  2.  Test  geometry  parameters. 


Parameter  Value 


Collector  1  Position  (Pc^) 
Collector  1  Velocity  (Vci) 
Collector  2  Position  (Pc2) 
Collector  2  Velocity  (Vc2) 
Emitter  Position  (Pg) 
Emitter  Velocity  (Vg) 


[0,  0,  7500]  m 
[100,  0,  0]  m/s 
[2000,  0,  7500]  m 
[100,  0,  0]  m/s 
[10000,  10000,  0]  m 
[0,  0,  0]  m/s 


Table  3.  Transmitted  signal  parameters. 


Parameter 

Value 

Carrier  Erequency  (/c) 

1000.025  MHz 

Sampling  Erequency  {fs) 

100  kHz 

Symbol  Rate  {Rsym) 

10000  symbol/sec 

Signal-to-Noise  Ratio  (SNR) 

Variable 

The  first  case  chosen  uses  a  signal  length  of  V  =  2^^  samples,  a  SNR  of  10  dB,  and  a  map 
grid  size  of  dm  =  100  meters.  Eurthermore  10  CAE  snapshots  were  computed  at  10  sec¬ 
ond  intervals  and  then  summed  to  provide  the  final  surface.  These  parameters  were  chosen 
because  the  computations  involved  are  short  enough  that  a  large  map  could  be  realistically 
generated  by  both  algorithms  and  a  large  scale  comparison  made.  Additionally,  this  case 
provides  a  reasonable  baseline  because  the  signal  length  is  small  enough  that  the  Hartwell 
algorithm  should  have  no  difficulty  and  will  perform  perfectly  fine.  As  shown  in  Eigure 
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16,  these  expeetations  are  realized.  Both  algorithms  peak  near  the  eorreet  loeation  of  the 
emitter  and  show  the  expeeted  ambiguity  about  the  y  axis  due  to  the  lead-trail  colleetor 
eonfiguration.  Significant  differences  can  be  noted,  however.  The  Moss  algorithm  provides 
a  much  smoother  surface  as  each  data  point  is  computed  independently  and  can  take  any 
value.  The  Hartwell  algorithm,  by  contrast,  maps  the  CAP  surface  to  the  geographic  plane 
by  placing  the  CAP  magnitude  values  into  xy  bins.  This  binning  behavior  is  not  necessar¬ 
ily  detrimental,  but  in  some  situations  it  could  lead  to  losing  some  valuable  information. 
Another  observation  of  note  from  this  test  is  that  the  smoothness  of  the  Moss  algorithm 
means  that  more  of  the  CAP  energy  is  located  away  from  the  peak.  Since  the  peak  is  what 
shows  the  emitter’s  location,  and  is  of  interest,  this  dispersion  could  make  geolocation  more 
difficult  with  a  greater  spread  in  the  Confidence  Interval  (Cl). 


Moss  algorithm  Hartwell  algorithm. 

Figure  16.  N  =  2^^,  dm  =  100,  SNR  =  10  dB,  10  Snapshots,  Full  Map. 

The  second  test  aims  to  see  the  effect  of  increasing  the  number  of  snapshots.  These 
snapshots  were  Hartwell’s  method  [5]  of  working  around  the  signal  length  limitations  of  his 
algorithm  by  spacing  out  collection  intervals  to  achieve  isodop  and  isochrone  rotation  with 
the  necessarily  short  collection  times.  In  principle,  this  method  of  dividing  the  collection 
into  snapshots  is  not  necessary  with  the  Moss  algorithm;  we  should  be  able  to  just  integrate 
straight  through.  However,  in  order  to  make  one-to-one  comparisons,  and  because  of  the 
computational  complexity  of  long  signals,  the  Moss  algorithm  will  also  use  snapshots  for 
these  tests.  Therefore,  this  test  leaves  all  parameters  equal  to  those  of  the  first  test,  but 
changes  the  number  of  snapshots  from  10  to  30. 
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As  shown  in  Figure  17,  the  results  are  not  entirely  as  expected.  For  the  Hartwell  al¬ 
gorithm,  the  peaks  are  indeed  much  narrower  as  desired,  but  the  map  also  becomes  much 
noisier.  This  is  attributed  two  factors.  First,  Hartwell’s  algorithm  only  computes  the  magni¬ 
tude  of  the  CAF,  which  is  then  non-coherently  summed,  losing  any  phase  information.  This 
leads  to  less  destructive  interference  occurring  and  the  noise  adding  up  with  each  snapshot. 
Second,  the  binning  characteristic  already  discussed  creates  an  ’on/off’  behavior  in  CAF 
values  at  given  xy  coordinates.  After  the  summation,  this  leads  to  some  xy  coordinates  hav¬ 
ing  very  large  values,  while  those  right  next  to  them  may  have  very  low  values.  The  Moss 
algorithm  does  not  show  this  characteristic.  The  surface  is  smooth,  without  any  binning, 
and  the  surfaces  are  all  summed  coherently  before  taking  the  surface  magnitude  as  the  final 
step  before  plotting.  This  seems  to  reduce  the  clutter  greatly  and  make  additional  snapshots 
more  beneficial.  This  implies  that  the  Moss  algorithm  would  perform  even  better  if  we  had 
the  computational  power  to  collect  straight  through,  without  snapshotting,  but  verifying 
this  is  left  for  future  work. 
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Figure  17.  N  =  2^^^,  dm  =  100,  SNR  =  10  dB,  30  Snapshots,  Full  Map. 


The  third  test  investigates  the  effect  of  increasing  the  SNR  to  30  dB,  and  also  doubles  the 
length  of  the  signal.  Additionally,  this  test  reduces  the  map  size  so  as  to  concentrate  on  just 
the  area  of  interest.  This  is  done  both  to  reduce  the  computation  time,  and  also  because  the 
first  two  tests  have  already  shown  the  large  scale  behavior,  which  mostly  follows  expecta¬ 
tions.  The  results  are  shown  in  Figure  18.  The  resolution  used  is  still  relatively  large,  and  at 
this  scale  both  algorithms  perform  adequately.  However,  it  can  be  seen  that  smoothness  of 
the  Moss  algorithm  again  makes  it  easier  to  find  the  true  peak  of  the  surface  and  determine 
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an  appropriate  Cl  for  the  geolocation  estimate.  The  Hartwell  algorithm  peaks  nearly  as 
close  to  the  true  emitter  location  as  the  Moss  algorithm,  but,  there  seems  to  be  a  greater 
variance  in  the  peak  and  multiple  locations  are  binned  into  the  same  value.  This  makes  it 
difficult  to  determine  the  true  peak. 
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Moss  algorithm  Hartwell  algorithm 

Figure  18.  N  =  2^^,  dm  =  100,  SNR  =  30  dB,  10  Snapshots. 


The  next  test  is  an  attempt  to  determine  if  the  Moss  algorithm  is  able  to  perform  better 
than  the  Hartwell  algorithm  for  long  duration  signals.  Achieving  this  was  the  primary 
motivation  for  this  work,  and  it  is  important  see  if  it  has  been  successful  or  not.  As  shown 
previously  in  Figure  11,  the  Hartwell  algorithm  begins  to  show  significant  smearing  and 
loss  of  accuracy  at  signal  lengths  of  2^^  so  a  length  of  2^^  was  chosen  for  this  test.  With 
signals  of  this  length  the  computation  time  becomes  a  serious  problem  so  very  small  map 
sizes  were  used  to  make  it  practical. 

The  results  of  this  test  are  shown  in  Figure  19.  The  first  thing  that  we  notice  is  that, 
this  time,  the  Moss  algorithm  is  just  as  rough  as  the  Hartwell  algorithm.  This  is  simply 
the  result  of  having  a  grid  resolution  on  the  same  order  as  the  grid  size,  the  next  test  will 
correct  this.  Beyond  the  roughness,  it  can  be  seen  is  that  both  algorithms  create  a  peak 
with  approximately  the  same  error.  It  appears  that  the  Moss  algorithm  has  a  slightly  lower 
(better)  variance,  but  this  is  impossible  to  determine  from  just  one  test  at  this  resolution,  so 
no  true  conclusions  can  be  drawn  from  this. 

In  order  to  draw  a  better  conclusion,  the  test  is  re-done  with  a  higher  resolution  of  dm  = 
10  meters.  This  increases  the  computation  time  by  a  factor  of  10,  so  the  grid  size  is  halved 
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Figure  19.  N  =  2^^,  dm  =  100,  SNR  =  30  dB,  10  Snapshots. 

to  bring  the  increased  computational  complexity  to  a  factor  of  4.  As  shown  in  Figure  20, 
the  results  appear  to  be  much  better.  The  Hartwell  algorithm’s  results  appear  smoothed  out 
some  and  it  arrives  at  a  closer  estimate  for  the  emitter’s  position.  However,  this  may  be 
deceiving.  The  surface  seems  to  be  monotonically  decreasing  away  from  the  peak,  but  this 
could  merely  be  a  result  of  using  such  a  small  grid  size,  and  the  true  peak  may  be  located 
off  the  grid.  With  this  resolution  and  grid  size  it  is  impossible  to  tell  if  this  is  the  true  peak, 
or  just  a  local  peak.  Furthermore,  even  if  it  is  the  true  peak,  the  energy  is  quite  dispersed 
over  the  grid,  which  will  result  in  a  higher  variance  in  the  maps  energy.  This  means  that  the 
statistical  geolocation  accuracy  will  most  likely  be  worse  if  multiple  runs  were  performed 
with  random  noise. 

The  Moss  algorithm,  on  the  other  hand,  is  able  to  show  the  true  oscillating  behavior 
of  the  CAF  surface,  which  gives  more  confidence  that  we  have  found  the  true  peak.  The 
energy  in  the  surface  is  also  more  concentrated  near  the  peak — meaning  that  if  more  tests 
could  be  completed,  and  statistical  analysis  performed,  the  results  are  more  likely  to  be 
consistent. 
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Moss  algorithm  Hartwell  algorithm 

Figure  20.  N  =  2^^,  dm  =  10,  SNR  =  30  dB,  10  Snapshots. 


4.3  Conclusion  of  Results 

Without  the  ability  to  perform  tests  with  both  large-seale  and  high  resolution  maps,  it  is 
diffieult  to  arrive  at  the  true  big-pieture  eharaeteristies  of  the  algorithm.  Additionally,  it  is 
impossible  to  conclusively  determine  the  accuracy  and  performance  without  the  ability  to 
do  multiple  runs  of  each  test  with  random  noise,  and  then  perform  statistical  analysis  to 
truly  see  what  the  mean  and  variance  of  the  detection  results  are.  However,  short  of  these 
conclusive  results,  it  appears  that  the  Moss  algorithm  succeeds  in  geolocating  RF  signals 
with  equivalent,  or  better,  accuracy  than  the  Hartwell  algorithm,  as  shown  in  Figures  16-20. 
Furthermore,  it  is  clear  that  the  Moss  algorithm  provides  a  much  better  mapping  resolution 
than  the  Hartwell  algorithm.  This  suggests  both  better  accuracy,  and  an  improved  ability  to 
distinguish  more  subtle  behaviors  that  would  otherwise  be  smoothed  out  and  missed.  More 
work  is  necessary  to  draw  convincing  conclusions  as  to  the  utility  of  the  Moss  algorithm, 
but  this  chapter  has  shown  that  it  merits  further  investigation.  The  final  chapter  will  discuss 
the  further  work  that  needs  to  be  accomplished  in  order  to  move  this  work  forward,  and  the 
current  thoughts  on  how  to  achieve  it. 
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CHAPTER  5: 
Future  Work 


This  work  has  demonstrated  that  the  Moss  CAFMAP  algorithm  is  capable  of  geolocating 
RF  emitters.  However,  there  is  still  a  lot  of  work  that  should  be  done  to  truly  complete  this 
study,  which  will  have  to  be  left  for  a  later  work  due  to  time  constraints.  This  section  will 
discuss  the  main  aspects  that  require  further  investigation. 

5.1  Reduce  Complexity 

In  order  to  be  practically  useful,  the  algorithm  must  be  able  to  run  in  less  time,  so  that 
results  can  be  obtained  in  near  real-time  as  required  by  real  world  operations.  The  first,  and 
simplest,  method  to  improve  the  performance  is  to  port  the  algorithm  to  a  more  suitable 
programming  language,  such  as  C  or  Fortran.  These  compiled  languages  are  much  better 
at  handling  large  computations  and  loops,  such  as  those  present  in  the  Moss  algorithm. 
MATLAB  was  used  in  this  work  due  to  its  ease  of  use,  and  the  availability  of  previous  code 
which  could  be  used  as  a  foundation.  However,  as  an  interpreted  language,  optimized  for 
linear  algebra,  it  does  not  perform  well  with  the  kinds  of  computations  required  in  this  algo¬ 
rithm.  It  is  expected  that  such  a  port  would  provide  substantial  performance  improvements, 
even  if  the  algorithm  itself  were  left  in  its  current  state. 

5.1.1  Adaptive  Resolution 

Adopting  an  adaptive  resolution  regime  is  another  method  to  reduce  the  complexity  of  the 
algorithm.  Rather  than  try  to  geolocate  the  emitter  in  a  single  run  it  could  be  done  in 
multiple  stages,  where  each  stage  adapts  the  grid  size,  grid  resolution,  and  signal  length  as 
appropriate  until  the  desired  accuracy  is  achieved.  The  initial  stage  could  use  a  very  short 
signal  length  and  large  grid  with  low  resolution  to  find  the  approximate  emitter  location, 
and  then  the  process  can  be  repeated  with  a  smaller  grid,  higher  resolution,  and  longer 
signal  length  to  iteratively  arrive  at  a  better  solution.  This  concept  is  very  similar  to  the 
coarse/fine  method  implemented  by  Johnson  [11]. 
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5.1.2  Parallelization 

One  great  advantages  of  modern  mieroproeessors  is  their  multi-eore  arehiteeture.  This  al¬ 
lows  them  to  simultaneously  perform  multiple  operations  to  reduee  proeessing  time.  The 
eurrent  CAFMAP  implementation  does  not  make  use  of  this  eapability,  and  is  a  faetor  that 
eould  be  greatly  improved.  The  algorithm  should  lend  itself  quite  well  to  parallelization,  as 
the  same  eomputations  are  performed  with  identieal  data  many  times  with  very  well  defined 
and  independent  ehanges.  This  has  been  shown  true  by  Singh  et  al.  as  they  sueeessfully  im¬ 
plemented  the  standard  CAF  in  the  49-eore  Maestro  proeessor,  with  signifieant  reduetions 
in  proeessing  time  [12].  It  should  be  fairly  simple  to  break  up  the  new  Moss  algorithm  into 
parallel  proeessing  streams  that  ean  be  reeombined  at  the  end.  For  example,  the  proeessing 
eould  be  broken  up  either  by  grid  loeations,  or  by  snapshots.  Both  ways  would  provide  a 
set  of  long  eomputations  that  are  entirely  independent  and  eould  be  performed  by  separate 
proeessors. 

5.1.3  Block  Updates 

In  the  eurrent  implementation  of  the  Moss  algorithm,  the  eolleetor-emitter  geometry  is  es¬ 
sentially  being  updated  with  every  sample.  While  this  refleets  the  true  reality,  it  is  useless 
in  praetiee.  At  orbital  veloeities  the  eolleetors  will  be  moving  approximately  7.5  km/see, 
and  with  a  representative  sampling  frequeney  of  100  kHz,  that  means  that  the  eolleetors 
will  ehange  position  by  km/sec^ =  0.075  m  between  eaeh  sample.  So  TDOA,  for 
example,  will  only  ehange  by  a  maximum  of  ^  ~  ~  ^  10^^®  see  between  sam¬ 

ples.  This  is  a  negligible  differenee,  that  is  rounded  out  by  the  algorithm  in  its  digital 
eomputation  anyway. 

Therefore,  another  approaeh  that  may  possibly  be  used  to  reduee  the  eomplexity,  would 
be  to  divide  the  eolleeted  signals  into  bloeks  of  eonstant  geometry,  and  only  update  eollee- 
tor  positions  when  moving  to  a  new  bloek.  By  itself  this  would  only  simplify  the  algorithm 
marginally,  but  it  may  make  it  possible  to  use  the  FFT  or  eorrelation  methods  to  ealeulate 
the  CAF  of  eaeh  bloek.  This  would  reduee  eomplexity  immensely,  as  those  methods  are 
far  more  effieient  than  brute  foree  ealeulations.  More  work  would  need  to  be  done  to  prove 
the  feasibility  of  using  the  FFT  or  eorrelation  funetions  however. 
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5.2  Experimentally  Determine  Statistical  Results 

Once  the  performanee  is  improved,  it  is  also  important  to  gather  statistieal  results  on  the 
aeeuraey  and  variations  in  the  algorithm.  It  is  impossible  to  eompletely  eharaeterize  and 
determine  the  usefulness  of  this  work  without  having  these  statisties  experimentally  deter¬ 
mined.  This  is  impossible  to  earry  out  with  the  eurrent  implementation,  due  to  the  eompu- 
tational  eomplexity  and  long  run  times,  but  onee  those  issues  are  resolved,  the  eharaeteri- 
zation  of  the  algorithm  needs  to  be  eompleted. 

5.3  Improve  Geometry 

Finally,  this  work  has  used  a  2-Dimensional  geometry  for  the  emitter,  with  linear  motion,  in 
order  to  keep  eomplexity  at  a  minimum  while  performing  the  proof  of  eoneept  for  the  work. 
However,  this  geometry  is  only  realistie  in  limited  situations,  sueh  as  aireraft  performing 
the  eolleetion  over  a  small  enough  land  area  that  ean  be  aeeurately  approximated  as  flat. 
For  more  general  applieability,  the  algorithm  must  be  fully  extended  to  three  dimensions, 
and  for  spaee  applieations,  the  eolleetors  must  follow  an  elliptieal  path  for  updating  the 
geometry.  Additionally,  it  is  desirable  to  allow  the  emitter  to  be  mobile,  as  well  as  the 
eolleetors,  sinee  this  is  often  the  ease  in  real-world  applieations. 

5.4  Conclusion 

In  eonelusion,  it  is  shown  that  the  reformulated  equations  for  TDOA  and  FDOA  are  ae- 
eurate,  and  ean  be  used  in  a  modified  CAFMAP  algorithm  to  aeeount  for  the  motion  of 
the  eolleetors  over  time.  The  resulting  Moss  algorithm  is  able  to  geoloeate  RF  emitters 
with  aeeuraey  similar  to  that  of  the  original  Hartwell  algorithm  for  shorter  eolleetion  times, 
and  with  apparent  greater  aeeuraey  for  long  eolleetion  times.  The  improvement  eomes, 
however,  with  the  eost  of  high  eomputational  eomplexity.  This  is  signifieant  enough  that  it 
must  be  mitigated  before  the  new  algorithm  ean  be  praetieally  implemented,  but  the  work 
required  is  merited.  Additionally,  the  approaeh  used  here,  of  modifying  the  TDOA  and 
FDOA  equations  to  aeeount  for  eolleetor  motion,  should  be  generally  applieable  to  many 
other  methods  of  geoloeation  that  estimate  TDOA  and  FDOA  as  a  fundamental  step.  This 
eould  lead  to  mueh  broader  use  of  these  results,  onee  the  eomputational  issues  have  been 
solved,  and  improve  the  geoloeation  aeeuraey  of  multiple  methods  in  eommon  use  today. 
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APPENDIX:  MATLAB  Code 


%  This  is  an  extension  of  Hartwell's  CAF_MAP  method  and  much  of  the 
%  initial  code  structure  was  taken  from  his  work. 

o, 

o 

%  This  function  defines  all  of  the  key  parameters  of  the  system, 

%  generates  the  signal  vectors  using  Joe  Johnson's  sig_gen  function, 
%  calls  caf_map  to  generate  the  CAF_MAP  in  X  and  Y  geographic 
%  coordinates,  and  then  plots  the  combined  result 

o, 

0 

%  Usage:  This  is  the  initiating  script,  no  inputs  or 
%  outputs  are  required. 

o, 

0 

%  Written  by:  LT  Andrew  Moss,  Jan.  2014 
%  Last  Modified:  LT  Andrew  Moss  December  2014 


%  clear  all;close  all;clc 


close 

all 

NUM_SNAPS  =30;  %  How 

GAP  = 

10;  %  Gap 

Pci  = 

[0,0,7.5e3];  % 

Vcl  = 

[100,0,0];  % 

Pc2  = 

[2e3,0,7.5e3];% 

Vc2  = 

[100,0,0];  % 

Pe  = 

[10000,10000, 0] ; 

Ve  = 

O 

o 

o 

fO  = 

1000.025e6;  % 

f  s  = 

100e3;  % 

Rsym 

=  10e3;  % 

N  =  2 

"21;  % 

%  Emitter  velocity  in  meters/sec. 


%  integration  time. 


Es_Nol  =  30; 


%  Received  SNR  for  collector  1. 
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Es_No2  =30;  %  Received  SNR  for  collector  2. 

dm  =  100;  %  Resolution  of  the  generated  XY  map  in  meteres. 

%  This  creates  the  XY  grid.  Starts  from  Pel  and  makes  a  square  to  Pe2 
%  with  a  resolution  of  dm  as  specified  above.  Units  are  meters. 

Pel  =  [  8 . 5e3, 8 . 5e3 ] ; 

Pe2  =  [11.5e3,11.5e3] ; 

t_snap  =  N/fs;  %  The  time  covered  by  each  snapshot  in  seconds. 

indexX  =  Pel  ( 1 )  : dm : Pe2 ( 1 ) ;  %  X  grid  points 
indexY  =  Pel  (2 )  : dm : Pe2 (2 ) ;  %  Y  grid  points 

map_comb  =  zeros (length (indexX) , length (indexY) ) ;  %  initialize  map. 

%  Break  up  the  signals  into  snapshots  and  find  the  caf_map  of  each  and 
%  sum  the  magnitudes 
for  ii  =  1 : 1 :NUM_SNAPS 

%  Update  Collector  1  position. 

Pcl_snap  =  Pci  +  ii*t_snap*Vcl  +  (ii— 1) *GAP*Vcl; 

%  Update  Collector  2  position. 

Pc2_snap  =  Pc2  +  ii*t_snap*Vc2  +  (ii— 1) *GAP*Vc2; 

%  Creates  the  signal  vectors  for  the  current  snapshot. 

[Sal,  Sa2,  Sl_snap,  S2_snap,  Peel,  Pcc2]  =  ... 
sig_gen (Pcl_snap, Vcl , Pc2_snap, Vc2 , . . . 

Pe, Ve, f 0 , f s , Rsym, N, Es_Nol , Es_No2 ) ; 

%  Call  the  caf_map  functions  to  generate  the  CAE  surface 
%  in  the  XY  plane 
%  for  the  current  snapshot. 

[map, indexX, indexY] =caf_map (Sl_snap, S2_snap, fO, . . . 

f s , dm, Pel , Pe2 , Pcl_snap, Vcl , . . . 

Pc2_snap, Vc2 )  ; 

%  Add  the  current  snapshot  to  the  sum  of  the  previous 
%  snapshots.  All  snapshots  will  have  a  peak  at  the  actual 
%  emitter  location  and  will  interfere  both  constructively 
%  and  destructively  elsewhere. 
map_comb  =  map_comb  +  map; 
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end 


%  Plotting  functions  that  may  be  used  if  found  useful. 

%  [tdoa_grid,  fdoa_grid,  PtempX, PtempY]  =  ... 

%  tdoa_fdoa_grid3d (Pcl_snap, Vcl , Pc2_snap, Vc2 ,Pel,Pe2,fO, dm) ; 

%  Tau_Lo  =  round (min (min (tdoa_grid) ) *fs) —10; 

%  Tau_Hi  =  round (max (max (tdoa_grid) ) *fs) +10; 

%  Freq_Lo  =  (min (min (fdoa_grid) ) /fs)— . 001; 

%  Freq_Hi  =  (max (max ( f doa_grid) ) /f s ) + . 001 ; 

%  [TDOA,  FDOA,  MaxAmb,  Amb,  TauValues , FreqValues ]  =  ... 

%  caf_peak ( Sl_snap,  S2_snap,  Tau_Lo, . . . 

%  Tau_Hi,  Freq_Lo,  Freq_Hi,  fs,  0)  ; 

%  figure ( ) 

%  f=mesh (TauValues/Fs, FreqValues*Fs, (abs (Amb) ) ) ; 

%  title ( 'CAF' ) ; 

%  set  ( f,  ' XData TauValues ) ; 

%  set  ( f,  ' YData FreqValues )  ; 

%  xlabel ( ' TDOA' ) ; 

%  ylabel ( 'FDOA' ) ; 

%  axis  tight 

savefilel  =  ' map_comb . mat ' ; 
save ( savef ilel ,  'map_comb'); 

%  Plot  mesh  of  the  combined  snapshots 
figure  ( ) 
hold  on 

f=mesh( (abs (map_comb' ) ) ) ; 
set  ( f ,  'XData ' , indexX) ; 
set  ( f ,  'YData ' , indexY) ; 
title ( ' ' ) 

xlabel ('X  (meters) ') 
ylabel ('Y  (meters) ') 
zlabelC  I  CAF  I  ') 
axis  tight 
grid  on 

set (gca, ' Font Size ' , 12 ) 

set ( f indall (gcf ,  ' type ' ,  ' text ' ) ,  ' Font Size  '  , 12 ) 


function  [CAF, indexX, indexY] = . . . 
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caf_map (Sl,S2,Fo,Fs, dm, Pel,Pe2,Pcl,Vcl,Pc2,Vc2) 

%-k'k-k-k-k-k-k-k-k-k'k-k'k'k-k-k-k-k-k-k-k'k-k'k-k-k'k-k-k-k-k-k'k-k'k'k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k 

%  [map, PtempX, PtempY] =caf_map (Sl,S2,Fo,Fs, dm, Pel,Pe2,Pcl,Vcl,Pc2,Vc2) 

%  This  function  will  calculate  a  CAF  surface  based  upon  input  signals 
%  SI  &  S2  and  map  the  caf  to  the  2  dimensional  plane  given  by 
%  Pel  and  Pe2 . 

%  Inputs; 

%  SI  Signal  from  collector  1. 

%  S2  Signal  from  collector  2  . 

%  Fo  Carrier  Frequency  of  signal. 

%  Fs  Sampling  Rate  . 

%  dm  Resolution  in  meters  . 

%  Pel  Start  of  grid  for  Emitter's  Position  [X,Y]  in  meters. 

%  Pe2  End  of  grid  for  Emitter's  Position  [X,Y]  in  meters. 

%  Pci  Position  of  collector  1  and  beginning  of  snapshot  [X,Y,Z] . 

%  Vcl  Velocity  Vector  of  collector  1  [x,y,z] . 

%  Pc2  Position  of  collector  2  and  beginning  of  snapshot  [X,Y,Z] . 

%  Vc2  Velocity  Vector  of  collector  2  [x,y,z] . 

9-9'9-9-9-S-9-9'9-9-9'9-9'9-9-9-9-9-9'9'9-9'9-9'9-9-9-2-9-9'9'9-9'9-9'9-2-9-9'9-9'9' 

oooooooooooooooooooooooooooooooooooooooooo 

%  Written  By:  LT  Andrew  Moss,  April  2014 
%  Last  Modified  By:  LT  Andrew  Moss,  December  2014 

%  'k-k-k'k-k-k-k-k-k'k-k'k'k-k-k-k-k-k-k-k'k-k'k-k-k'k-k-k-k-k-k'k-k'k'k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k'k-k-k'k-k-k-k-k'k-k-k'k-k-k-k-k 

c  =  2.997925e8;  %  Speed  of  light  in  m/s . 

Ts  =  1/Fs;  %  Sampling  period  in  sec. 

N  =  length (SI )  ; 

%  Builds  the  position  vectors  for  the  Emitter's  Position 
%  Note  this  assumes  a  flat  earth  and  the  emitter  is  a  Om  alt. 
indexX  =  Pel(l) :dm:Pe2(l);  %  X  grid  points  in  meters 

indexY  =  Pel  (2)  :dm:Pe2(2);  %  Y  grid  points  in  meters 

%  Get  the  dimensions  of  the  map  in  order  to  go  through  each 

%  location. 

Nx  =  length ( indexX) ; 

Ny  =  length ( indexY) ; 

%  Initialize  matrices  for  speed. 

CAF  =  zeros (Nx, Ny) ; 

TauValue  =  zeros (1,N); 
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FDOA 


zeros ( 1 , N) ; 


%  Start  going  through  each  map  location  and  calculated  the  CAF  value  if 
%  that  grid  were  the  actual  location  of  the  emitter, 
for  11=1: Nx 

for  j  j  =  1 : Ny 

%  Directly  computing  A(x,y)  =  sum{  si (t ) s2 * (t , x, y ) exp  {  j2pi  t 
%  f(t,x,y))  as  derived  in  notes  on  160CT13 
m  =  zeros { 1 , N) ; 
exponent  =  zeros (1,N); 

Pe  =  [ indexX ( ii ) , indexY ( j j ) , 0 ] ; 

S2_conj  =  conj(S2); 

%  Loop  through  each  element  of  S2  for  a  given  map  location, 
for  nn  =  l:length(S2) 

%  TauValue  is  calculated  as  a  function  of  emitter  position, 

%  collector  geometry  and  collector  velocity.  Must  be 
%  rounded  to  a  whole  number. 

%  Calculate  relative  position  vectors  for  each  collector  and 
%  the  current  map  location.  Update  each  sample  to  account 
%  for  the  position  change  of  the  collectors  over  time, 
rl  =  Pe  —  (Pcl+Vcl . *nn*Ts ) ; 
r2  =  Pe  —  (Pc2+Vc2 . *nn*Ts )  ; 

%  Find  TDOA. 

TauValue (nn)  =  round ( (norm { r2) —norm (rl ))/ (c*Ts ))  ; 

%  Perform  the  shift  in  the  collected  signals  to  account 
%  for  the  TDOA  and  realign  the  signals.  Then  calculate 
%  the  mixing  function. 

if  TauValue (nn)  >=  0  &&  (nn+TauValue (nn) )  <  N 
m(nn)  =  SI  (nn)  *  S2_con j (nn+TauValue (nn) )  ; 
elseif  (nn— TauValue (nn) )  <  N 

m(nn)  =  SI (nn— TauValue (nn) )  *  S2_conj (nn) ; 

end 

%  Calculate  the  FDOA  from  the  relative  position  vectors. 

%  This  are  done  as  vectors  to  avoid  more  loops,  but  as  with 


43 


%  TDOA  the  value  is  changing  in  time. 

FDOA(nn)  =  Fo/c* {Vc2 *r2 ' /norm (r2) —Vcl*rl ' /norm (rl )) ; 

%  Calculate  the  exponential  term  of  the  CAF . 
exponent (nn)  =  exp (— 1 j*2*pi*FD0A (nn) /Fs* (nn— 1) )  ; 

end 

%  Finish  calculating  the  CAF  value  for  the  current  map  position, 
product  =  m. *exponent; 

CAF(ii,jj)  =  sum (product ) ; 


%  The  following 


f doa_grid2 ( ii, j j )  =  mean(FDOA); 
tdoa_grid2 (ii, j j )  =  mean (TauValue) ; 

T_vals  =  linspace (min (min (tdoa_grid2)  )  ,  .  .  . 

max (max (tdoa_grid2 ) ) , length 

%%  Just  Calc  the  CAF  and  then  map  it  via  Hartwell's  method 
%  adds  the  3rd  dimensions  at  0  meters  in  altitude 
gridP  =  [ indexX (ii ) , indexY ( j j )  ,  0 ]  ; 


dopplerl ( ii , j j )  = 


doppler2 ( ii , j j )  = 


Fo/c  *  dot(Ve— Vcl,  gridP— Pci) 
/  norm(gridP  —  Pci); 

Fo/c  *  dot(Ve— Vc2,  gridP— Pc2) 
/  norm(gridP  —  Pc2); 


%  Calculates  the  FDOA 


fdoa_grid (ii, j j )  =  dopplerl ( ii , j j )  —  doppler2 ( ii , j j ) ; 


%  Calculates  the  TDOA 

tdoa_grid (ii, j j )  =  — (norm(gridP  —  Pc2)  — ... 

norm(gridP  —  Pci))  /  c; 


end 


end 
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save file 1  =  'CAF.mat'; 
save (savefilel,  ' CAF ' ) ; 


function  [Sal,  Sa2,  SI,  S2,  Peel,  Pcc2]  =  sig_gen (Pel , Vcl , Pc2 , Vc2 , . . . 

Pe, Ve, f 0 ,  f s , Rsym, N, Es_Nol , Es_No2 ) 

%  -k-k-k'k-k-k-k-k-k'k-k'k'k-k'k-k-k-k-k-k-k-k-k'k'k'k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k'k-k-k-k-k-k'k-k 

%  [Sal,  Sa2,  SI,  S2,  Pel,  Pc2]]  =  sig_gen; 

%  S1G_GEN  generates  BPSK  signal  pairs  based  upon  user— defined  param— 

%  eters  and  Cartesian  emitter— collector  geometries.  There  ari 

%  no  input  arguments,  since  the  function  queries  the  user  for 

%  all  required  inputs.  The  function  returns  four  vectors: 

%  Sal,  Sa2,  SI  &  S2.  These  are  the  Analytic  Signal  represen— 

%  tations  of  the  two  generated  signals,  and  the  Real  represen- 

%  tations  of  the  two  signals,  respectively. 

o, 

o 

%  Written  by:  LCDR  Joe  J.  Johnson,  USN 
%  Last  modified:  28  June  2005  By  Glenn  D.  Hartwell 
%  Modified  for  3D  simulations  and  export  collectors  positions 

%  -k'k-k'k-k-k-k-k-k'k-k'k'k-k'k-k-k-k-k-k-k-k-k'k'k'k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k-k-k-k-k'k-k-k-k-k-k-k'k-k-k'k'k-k-k-k-k-k-k-k-k-k-k'k-k 


clc; 

%  disp ( '  ' ) ; 


%  disp {'All  positions  and  velocites  must  be  entered  in  vector  format,'); 
%  disp('e.g.,  [X  Y  Z]  or  [X,  Y,  Z]  (including  the  brackets).'); 

%  disp { '  ' ) ; 

o, 

o 

%  Pcl(l,:)  =  input... 

%  ('Collector  l''s  POSITION  Vector  at  time  0  (in  meters)?  '); 

%  Vcl  =  input (' Collector  l''s  VELOCITY  Vector  (in  m/s)?  ');  disp('  '); 

o, 

o 


%  Pc2(l,:)  =  input... 

%  ('Collector  2''s  POSITION  Vector  at  time  0  (in  meters)?  ' ) ; 

%  Vc2  =  input (' Collector  2''s  VELOCITY  Vector  (in  m/s)?  ');  disp('  '); 

o 

0 


%  Pe  ( 1 ,  : )  =  input .  . . 

%  ('Emitter' 's  POSITION  Vector  at  time  0  (in  meters)?  '); 
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%  Ve  =  input (' Emitter '' s  VELOCITY  Vector  (in  m/s)?  ');  disp{'  '); 

o, 

0 

%  %  fO  and  fs  are  the  same  for  BOTH  collectors! 

%  fO  =  input (' Carrier  Frequency  (in  Hz)?  '); 

%  fs  =  input (' Sampling  Frequency  (in  Hz)?  '); 

Ts  =  1/fs;  %  Calculates  Sample  Period 

o, 

o 

%  Rsym  =  input (' Symbol  Rate  (in  symbols/s)?  ' ) ;  disp('  '); 

Tsym  =  1/Rsym;  %  Calculates  Symbol  Period 

o, 

0 


%  N  =  input ('How  many  samples?  ');  disp ( '  '); 

o, 

0 


%  Es_Nol 

=  input (' Desired 

Es/No 

at 

Collector  1 

(in  dB) ?  '  )  ; 

%  Es_Nol 

o. 

=  10^ (Es_Nol/10) ; 

%  Converts  from  dB  to  ratio 

0 

%  Es_No2 

=  input (' Desired 

Es/No 

at 

Collector  2 

(in  dB) ?  '  )  ; 

%  disp ( ' 

' ) ; 

Es_No2  = 

10" (Es_No2/10) ; 

%  Converts 

from  dB  to  ratio 

Pci 

=  [Pci;  z€ 

^ros 

(N-1, 

3)  ]  ; 

%  Initializing  all  the  matrices  makes 

Pel 

=  zeros (N, 

3) 

r 

%  later  computations  much  faster. 

Pc2 

=  [Pc2;  z€ 

^ros 

(N-l, 

3)  ]  ; 

Pe2 

=  zeros (N, 

3) 

r 

tl  = 

zeros  ( 1 , 

N)  ; 

t2  = 

zeros  ( 1 , 

N)  ; 

SI  = 

zeros ( 1 , 

N)  ; 

S2  = 

zeros  ( 1 , 

N)  ; 

A  =  1;  %  Amplitude  of  Signal 

c  =  2.997925e8;  %  Speed  of  light  in  m/s 
Ps  =  (A^2)/2;  %  Power  of  Signal 

sigmal  =  sqrt (Ps*Tsym/Es_Nol ) ;  %  Calculate  Noise  Amplification  fac— 

sigma2  =  sqrt (Ps*Tsym/Es_No2 ) ;  %  tors  using  Es/No  =  Ps*Tsym/sigma^2 


Noisel  =  sigmal . *randn (N,  1); 

Noise2  =  sigma2 . *randn (N,  1); 


%  Random  Noise  values  for  Signal  1 
%  Random  Noise  values  for  Signal  2 
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%  Builds  the  position  vectors  for  the  two  collectors 
for  index  =  2  :  N 

Pci { index, : )  =  Pci ( index  —  1 , : )  +  Ts*Vcl ; 

Pc2 { index, : )  =  Pc2 ( index  —  1 , : )  +  Ts*Vc2 ; 

end 


%  While  loop  determines  first  elements  of  Pel  and  tl.  tl(l)  is  the 
%  time  AT  THE  EMITTER  that  produces  the  1st  sample  received  at 
%  collector  1!  Pel(l,:)  is  the  position  of  the  emitter  when  it 
%  produces  the  1st  sample  received  by  collector  1. 

temp  =  inf;  %  Ensures  while  loop  executes  at  least  once 

tl(l)  =  0; 
tempPe  =  Pe ( 1 ,  : )  ; 
while  abs (temp  —  tl(l))  >  1/fO 
temp  =  tl  ( 1 ) ; 

tl(l)  =  —norm (Pci ( 1 ,: )  —  tempPe)  /  c; 
tempPe  =  Pe(l,:)  +  tl(l)*Ve; 

end 

Pel(l,:)  =  tempPe; 


%  While  loop  determines  first  elements  of  Pe2  and  t2.  t2(l)  is  the 

%  time  AT  THE  EMITTER  that  produces  the  1st  sample  received  at 
%  collector  2!  Pe2(l,:)  is  the  position  of  the  emitter  when  it 

%  produces  the  1st  sample  received  by  collector  2. 

temp  =  inf;  %  Ensures  while  loop  executes  at  least  once 

t2(l)  =  0; 
tempPe  =  Pe  ( 1 ,  :  )  ; 
while  abs (temp  —  t2(l))  >  1/fO 
temp  =  t2  ( 1 ) ; 

t2(l)  =  — norm (Pc2 ( 1 , : )  —  tempPe)  /  c; 
tempPe  =  Pe(l,:)  +  t2(l)*Ve; 

end 

Pe2(l,:)  =  tempPe; 
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%  Platform  positions  at  middle  of  snapshot 
Pccl= (Pci (N/2,  : ) )  ; 

Pcc2= (Pc2 (N/2, : ) ) ; 

%  Determines  the  earliest  time  at  the  emitter  for  this  pair  of  signals. 
StartPoint  =  min(tl(l),  t2(l)); 


%  Next  2  lines  determine  offsets  needed  for  signals  1  &  2  to  enter  the 
%  phase  vector  (P) .  This  simply  ensures  proper  line  up  so  that  bit 
%  changes  occur  at  the  right  times. 

Symbollndexl  =  1  +  f loor (abs (tl ( 1 )  —  t2{l))/Tsym)  *  (tl ( 1 ) >t2 { 1 ) ) ; 
SyinbolIndex2  =  1  +  floor  (abs  (tl  ( 1 )  —  t2(l))/Tsym)  *  (t2  ( 1 )  >tl  ( 1 )  )  ; 

for  index  =  2  :  N  %  Builds  the  Pel  and  tl  vectors 

temp  =  inf; 
tl (index)  =  0 ; 

%  1st  guess  is  that  emitter  will  advance  exactly  Ts  seconds. 
tempPe  =  Pel  (1,  : )  +  (tl (index  —1)  +  Ts) *Ve; 

%  While  loop  iteratively  determines  actual  time  &  position  for 
%  emitter,  based  on  instantaneous  geometry, 
while  abs (temp  —  tl (index) )  >  1/fO 

temp  =  tl (index); 

tl (index)  =  (index  —  l)*Ts  —  norm (Pci ( index, : )  —  tempPe)  /  c; 

%  Due  to  negative  times,  must  multiply  Ve  by  ELAPSED  time! 
tempPe  =  Pel(l,:)  +  abs (tl (1) —tl (index) ) *Ve; 
end 

Pel (index, :)  =  tempPe; 

end 


for  index  =  2  :  N  %Builds  the  Pe2  and  t2  vectors 

temp  =  inf; 
t2 (index)  =  0 ; 

%  1st  guess  is  that  emitter  will  advance  exactly  Ts  seconds. 
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tempPe  =  Pe2  (1,  :  )  +  (t2 (index  —1)  +  Ts) *Ve; 


%  While  loop  iteratively  determines  actual  time  &  position  for 
%  emitter,  based  on  instantaneous  geometry, 
while  abs (temp  —  t2 (index) )  >  1/fO 

temp  =  t2  (index) ; 

t2 (index)  =  (index  —  l)*Ts  —  norm (Pc2 ( index,  : )  —  tempPe)  /  c; 

%  Due  to  negative  times,  must  multiply  Ve  by  ELAPSED  time! 
tempPe  =  Pe2(l,:)  +  abs (t2 (1) —t2 (index) ) *Ve; 

end 

Pe2 (index, :)  =  tempPe; 

end 


%  Could  change  this  seed  to  whatever  you  want;  or  could  have  user 
%  define  it  as  an  input.  This  just  ensures,  for  simulation  purposes 
%  that  every  time  the  program  is  run,  the  BPSK  signals  created  will 
%  have  the  same  random  set  of  data  bits, 
rand ( ' seed ' ,  5 )  ; 

%  Create  enough  random  #'s  to  figure  phase  shift  (data  bits) 
r  =  rand (N, 1 ) ; 

P  =  (r  >  0.5)*0  +  (r  <=  0.5) *1;  %  Since  BPSK,  random  #  determines 

%  if  phase  is  0  or  pi 


%  Building  Xmitted  Signal  #1  vector...  These  represent  the  pieces  of 
%  the  signal  that  were  transmitted  by  the  emitter  to  arrive  at 
%  Collector  1  at  its  sample  intervals. 

Sl(l)  =  A*cos (2*pi*f 0*tl (1)  +  P ( Symbollndexl ) *pi )  +  Noisel(l); 

%  The  if  statement  inside  the  loop  changes  the  data  bit  if  the  time 
%  has  advanced  into  the  next  symbol  period, 
for  index  =  2  :  N 

if  tl (index)  —  StartPoint  >=  (Symbollndexl )  *  Tsym 

Symbollndexl  =  Symbollndexl  +  1; 

end 
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SI  (index) 

end 


A*cos (2*pi*f0*tl (index) 
Noisel (index) ; 


P ( Symbollndexl ) *pi )  + 


Sal  =  hilbert (SI)  ;  %  Calculates  the  ANALYTIC  SIGNAL  of  SI.  To 

%  compute  the  COMPLEX  ENVELOPE,  multiply  Sal 
%  by  . *exp (— j *2*pi*f 0 . *tl ) ; 


%  Building  Xmitted  Signal  #2  vector...  These  represent  the  pieces  of 
%  the  signal  that  were  transmitted  by  the  emitter  to  arrive  at 
%  Collector  2  at  its  sample  intervals. 

S2(l)  =  A*cos (2*pi*f 0*t2 (1)  +  P ( SymbolIndex2 ) *pi )  +  Noise2(l); 

%  The  if  statement  inside  the  loop  changes  the  data  bit  if  the  time 
%  has  advanced  into  the  next  symbol  period, 
for  index  =  2  :  N 

if  t2 (index)  —  StartPoint  >=  ( Symbol Index2 )  *  Tsym 

SymbolIndex2  =  SymbolIndex2  +  1; 

end 

S2 (index)  =  A*cos (2*pi*f0*t2 (index)  +  P ( SymbolIndex2 ) *pi )  +  ... 

Noise2 (index) ; 

end 

Sa2  =  hilbert (S2);  %  Calculates  the  ANALYTIC  SIGNAL  of  S2 .  To 

%  compute  the  COMPLEX  ENVELOPE,  multiply  Sa2 
%  by  . *exp (— j *2*pi*f 0 . *t2 ) ; 


%  This  function  call  simply  calculates  and  displays  the  expected  TDOAs 
%  and  FDOAs  at  the  Beginning  and  End  of  the  collection  time. 


%  tdoa_fdoa (fO, Pel (1, : ) ,Pel (N, : ) ,Pe2 (1, : ) ,Pe2 (N, :), Ve, Pci ( 1 ,:),.. . 
%  Pci (N, : ) , Vcl, Pc2 (1, ; ) , Pc2 (N, : ) , Vc 
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